Reducing eccentricity in black-hole binary evolutions with initial parameters from 

post-Newtonian inspiral 



Sascha Husa, Mark Hannam, Jose A. Gonzalez, Ulrich Sperhake, Bernd Briigmann 

Theoretical Physics Institute, University of Jena, 07743 Jena, Germany 

Standard choices of quasi-circular orbit parameters for black-hole binary evolutions result in 
eccentric inspiral. We introduce a conceptually simple method, which is to integrate the post- 
Newtonian equations of motion through hundreds of orbits, and read off the values of the momenta 
at the separation at which we wish to start a fully general relativistic numerical evolution. For the 
particular case of non-spinning equal-mass inspiral with an initial coordinate separation of D — 11M 
we show that this approach reduces the eccentricity by at least a factor of five to e < 0.002 as 
compared to using standard quasi-circular initial parameters. 
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I. INTRODUCTION 

Recent breakthroughs in numerical relativity P, 0, 0] 
have made it possible to accurately simulate the last or- 
bits, merger and ringdown of a black-hole binary sys- 
tem, and to compute the gravitational waves emitted in 
the process. Comparison of these waveforms with those 
produced by analytic techniques (i.e., post-Newtonian 
[PN] methods) has already begun 0, [!, 0, 0, i, Q, as 
has the process of preparingthese waveforms for use in 
gravitational- wave searches IS;, 9]. We would like to start 
the comparison and template-building process with bi- 
naries that undergo non-eccentric inspiral, but many of 
these simulations model systems with non-zero eccentric- 
ity. This is ironically due to the use of "quasicircular 
orbit" initial momenta, which are based on an approxi- 
mate circularity condition. However, a spiraling motion 
must contain a radial component. We wish to find initial 
tangential and radial momenta that, for a given initial 
separation, lead to non-eccentric inspiral. 

Pfciffcr, et al. recently suggested an iterative pro- 
cedure to reduce eccentricity. They simulate a system 
with quasi-circular (QC) parameters for two orbits, mea- 
sure the eccentricity, make an appropriate modification 
to the initial parameters (including the introduction of a 
radial component to the motion), and start the simula- 
tion again. They repeat this procedure until the eccen- 
tricity is reduced by a factor of ten. The drawback of 
their method is that it requires at least one "false start" , 
which is computationally expensive, and, as pointed out 
in [lfj, it is not clear how to generalize the method to 
evolutions of spinning black holes, for which the black- 
hole separation and wave frequency will not in general 
be a monotonic function of time. We would rather find a 
general method to calculate low-eccentricity parameters 
from the outset. 

In earlier papers we used a PN approximation to cal- 
culate QC parameters for equal-mass 11], unequal-mass 
[l2| and spinning [l3| black holes, and found that they 
compared well with parameters calculated using more so- 
phisticated numerical methods [ill ]. In this paper we 
follow up on Miller's work [l4| and use PN methods to 



estimate initial parameters for low-eccentricity inspiral. 

Our procedure is to numerically integrate the PN equa- 
tions of motion (at the highest PN order available; see for 
example fl5j ) for two point particles over hundreds of or- 
bits, and read off the particles' momenta when they have 
reached the separation we wish to use as an initial sepa- 
ration in a fully general relativistic numerical simulation. 
We use Mathematica for the integration, which typically 
takes several seconds. Note that we could instead use 
lower-order PN estimates of the radial momentum (for 
example by using the quadrupole formula), but we wish 
to obtain the most accurate results possible, and to use a 
method that can later be applied with some confidence to 
spinning binaries. A PN approach was also used to intro- 
duce a radial component to the motion in [l6|, although 
the details of the method were not given. 

The key to this approach is to exploit the fact that any 
initial eccentricity will decay over time due to the circu- 
larizing effect of gravitational-wave emission, but on a 
time-scale of hundreds of orbits, not the < 10 orbits typ- 
ically simulated by a numerical code. We therefore start 
the PN evolution at sufficiently large initial separation 
D (in practice D = 40M) to allow radiation-reaction to 
circularize the orbits. (We will quote time and length in 
units of the total initial black-hole mass M; see [llj.) 

Given the parameters from PN inspiral, full GR nu- 
merical evolutions are then performed with the BAM 
code [ill E3 > usm g the moving-puncture method @, 3 to 
evolve Bowen-York initial data [18| in puncture form [19| 
and generated by a pseudo-spectral method [201 ] . In or- 
der to perform long evolutions with sufficient accuracy for 
the present purpose, it was crucial to modify our previous 
evolution algorithm to use sixth order accurate derivative 
operators [211 ] , instead of the more standard fourth-order 
accurate choice. 

We find that the PN-inspired initial momenta lead to 
evolutions with at least five times less eccentricity than 
their QC counterparts. 

In Section [TT] we summarize the PN equations that we 
use and the method to integrate them. In Section HTT1 we 
present results from simulations of an equal-mass binary. 
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II. INTEGRATION OF THE PN EQUATIONS 
OF MOTION 



We have used the PN equations of motion as described 
in [l5[ in the ADMTT gauge. We have implemented both 
the usual Taylor-expanded and the effective-one-body 
(EOB) versions of the Hamiltonian; the calculations pre- 
sented here are however all based on the Taylor-expanded 
version. The PN solution in the ADMTT gauge for a 
two-body system agrees with our Bowen-York puncture 
initial data up to 2PN order (see, for example, the explicit 
solutions in Appendix A of [22]). The conservative part 
of the Hamiltonian is given up to third PN or der, and 
was originally derived in [13, [H, [24[ , see also (2f| [2(| [27| • 
Radiation-reaction flux terms are calculated up to 3.5PN 
order beyond the quadrupole order, which is achieved 
by averaging the radiation flux over one orbit, assuming 
quasi-circular inspiral [H, [H, 30] . We have also included 
the leading-order spin-spin and spin-orbit coupling terms 
for the conservative part of the Hamiltonian [3ll . |32| . 1331 . 
and spin-induced radiation flux terms as described in [151 ] 
(and again averaged over one orbit). 

In the nonspinning case the PN equations of motion are 
a system of six coupled ordinary differential equations of 
the form 



dx l 
~~dt 

dpi 
dt 



dH 

dpi ' 
dH 

dx l 



F, 



(1) 
(2) 



where H is the PN-Hamiltonian (responsible for the con- 
servative part of the dynamics), x l is the separation vec- 
tor between the two particles and p 1 is the momentum of 
one particle in the ccnter-of-mass frame. In the spinning 
case the system is augmented by the evolution equations 
for the spins. The quantity Fi is the radiation-reaction 
flux term. We have used the evolution equations precisely 
in the form presented in (lEI ] . 

Starting at a suitably large initial separation (D = 
40 M is used in practice for an equal-mass binary), initial 
momenta are chosen using the 3PN formula given in 
We have checked that D — 2QM would be somewhat too 
close — small oscillations in the radius are still visible at 
D = UM, where we wish to start the evolution of the 
full Einstein equations; see Fig. (fT]). Integrating the PN 
equations from D = lOOAf (2071.5 orbits) makes very 
little difference for the inspiral parameters. For the full 
numerical evolution to start at D = UM the tangential 
component of the momentum would change by 4x 10 %, 
the radial component by 0.2% as compared to using a 
PN-inspiral from D = 40M . 

The equations are then integrated in Mathematica 5.2 
using the NDSolve function with different options for the 
integration algorithm, tolerance levels and internal pre- 
cision to check the accuracy of the results. Mathematica 
stops the integration automatically when the PN equa- 
tions of motion become ill-defined. Several consistency 
checks are applied to make sure the correct equations are 
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FIG. 1: The radial momentum component is plotted ver- 
sus separation for PN-inspirals starting from D = 20M and 
D = 40M. A separation of D = 20M is clearly not sufficient 
to produce non-eccentric inspiral parameters, since small os- 
cillations can still be seen at D — UM, while for D = 40M 
the initial eccentricity has essentially decayed away. 



integrated: When radiation reaction is switched off, en- 
ergy and angular momentum should be constant up to 
numerical error. To check this, we integrate the conser- 
vative equations of motion with initial separations 40M 
and 50M for about 200 orbits and monitor the relative 
decay of energy and angular momentum, which remain 
below 3.5 x 10~ 6 for the NDSolve options we have used 
for the results presented here. When radiation reaction 
is switched on, the PN equations of motion imply the 
identity 



dE _ dip IdJ 
~dt ~ ~dt ' \~d~t 



in the circular equal- mass case [15j, where ip is the or- 
bital phase, A the time independent unit vector in the 
direction of the orbital angular momentum, and (•) de- 
notes the orbital average. We have checked that for the 
actual inspiral (which is not exactly circular) this equa- 
tion (without taking the orbital average) is satisfied to 
better than 2 x 10~ 4 over the entire inspiral. 

The system displays some initial eccentricity, but this 
decays and by the time the particles are at a separation 
suitable for a numerical evolution (i.e., D < 20M), the 
inspiral has negligible eccentricity, as shown in Fig. (fT]). 
Similar plots are also shown in [14|. 

We now wish to perform a full GR numerical evolution 
of the last orbits of the binary system. The puncture ini- 
tial data solver requires as input the black hole's masses, 
positions, and momenta. Given the masses and some de- 
sired initial separation, we can read off the appropriate 
momenta from the integrated solution (x l (t) , p % (t)) of the 
PN equations of motion. 
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Configuration 


P x /M 


Py/M 


ec 




QC11 


T0.0899395 





0.012 


0.01 


Ell 


T0.0900993 


^7.09412 x 1CT 4 


0.002 


0.002 



TABLE I: Initial physical parameters for a standard "quasi- 
circular orbit" (QC11) and PN-inspired low-eccentricity (Ell) 
configuration. Both have an initial coordinate separation of 
D — 11M and the punctures are placed at y — ±5.5M. The 
initial eccentricities eo and are estimated using Eqns. © 
and (01). 



III. NUMERICAL RESULTS 

We consider the configurations shown in Table HI The 
black-hole punctures are placed at an initial separation 
of D = 11M. They are given either quasi-circular orbit 
(QC) parameters as estimated from the 2PN-accurate 
expression in [34], [35[ , or the PN-inspired low-eccentricity 
parameters described in Section [TTJ When evolved both 
initial configurations lead to about seven orbits before 
merger. 

In the notation of [ll|, the data were evolved with a 
grid setup of x?7=2 [5 x 64 : 5 x 128 : 6] using the sixth-order 
accurate spatial finite differencing stencils, as described 
in detail in [211 ] . Lower-resolution runs and convergence 
tests show that the simulations are cleanly sixth-order 
convergent up to around lOOOAf of evolution time and 
drop slightly in convergence order after that. In this 
paper we only need the simulation up to t = 1000M, at 
which time the uncertainty in D(t) is 0.6%; for all earlier 
times it is lower. 

Figure [2] shows the coordinate separation of the punc- 
tures as a function of time, for simulations with the 
QC11 and Ell initial parameters. The figure begins at 
t = 257M, the time at which the binary completes one 
orbit. Before that time there are oscillations in D(t) due 
to gauge adjustments; the initially stationary punctures 
pick up speed, the lapse and shift adapt to the dynamical 
gauge conditions, and the numerical grid points rapidly 
retract from the extra asymptotically flat ends in the 
puncture initial data [111 [33] • All of these effects preclude 
a meaningful estimate of the eccentricity during the first 
orbit. The figure ends at t = 1050M, when the system 
has completed a further four orbits. We can clearly see 
oscillations due to eccentricity in the QC11 data, while 
the Ell data appears relatively eccentricity-free. 

We use two methods to estimate the eccentricity, as 
also used in 

HUES- Assume that we know the zero- 
eccentricity quasi-circular inspiral for our system, and 
denote the corresponding coordinate separation of the 
punctures as a function of time as D c (t) and the orbital 
frequency as uj c (t) . The coordinate separation and orbital 
frequency for any given numerical evolution are D(t) and 
Lo{t). The eccentricity can be estimated by extrema in 



either 



or 



e D (t) 



D(t) - D c (t) 
D c (t) ■ 

Uj(t) - LO c (t) 



(3) 



(4) 



In practice we estimate D c (t) by fitting a curve through 
the numerical D(t) for the low-eccentricity simulation 



D c (t) = aT 1/2 + bT + cT 3/2 + dT 2 



(5) 



where T = Tm — t and Tm is a rough estimate of the 
merger time. For the Ell simulation, we choose tia = 
1270M. Similarly we follow 5] and fit a fourth-order 
polynomial in time through the u>(t) curve for the low- 
eccentricity simulation and obtain u> c (t). 
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FIG. 2: Coordinate separation of the punctures as a function 
of time for the quasi-circular (QCU) and PN-inspired low- 
eccentricity (Ell) initial parameters. 



For the quasi-circular simulation QC11, Eqn. ([3]) gives 
e D = 0.012 ±0.002 and Eqn. Q gives e w = 0.01 ±0.001. 
The uncertainties are estimated by repeating the calcula- 
tion with curves D c (t) and u) c (t) fit through the eccentric 
QC11 data. Note that the frequency method (|4|) gives a 
lower value than the separation method ((3|); similar re- 
sults were found in 0, [H| • 

For the low-eccentricity Ell simulation, we find eo 
(1.002 ± 0.001 and e u = 0.002 ± 0.0005. The large un- 
certainty in the value from the separation method is due 
to the larger uncertainty in the curve fit through D(t). 
For both estimates, however, the one firm conclusion we 
can draw is that the eccentricity in the Ell simulation 
is significantly lower (by a factor of at least five or six) 
than that for the QC11 simulation. 
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The functions eo{t) and e w (f) are shown in Figure [3] 
For the QC11 simulations we see clear oscillations due 
to the eccentricity. The curves for the Ell run are much 
noisier. This may be due to errors in the curve fit through 
the Ell D(t) and ui(t) being of a similar magnitude to 
the oscillations due to the remaining eccentricity. 
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FIG. 3: The functions e D (t) and e u (t) for the QC11 and Ell 
simulations. The extrema of these functions give an estimate 
of the eccentricity, as described in the text. 



separation/M 


-r (xlO -3 ) 


-Pr/M (xlO -3 ) 


Pt/M 


8.0 


5.3857 


2.0906 


0.112349 


8.5 


4.4944 


1.7023 


0.107614 


9.0 


3.7839 


1.4019 


0.103376 


9.5 


3.2133 


1.1670 


0.099561 


10.0 


2.7512 


0.9813 


0.096109 


10.5 


2.3736 


0.8328 


0.092968 


11.0 


2.0624 


0.7128 


0.090099 


11.5 


1.8039 


0.6150 


0.087464 


12.0 


1.5872 


0.5343 


0.085037 


12.5 


1.4042 


0.4672 


0.082791 


13.0 


1.2487 


0.4110 


0.080706 


13.5 


1.1156 


0.3635 


0.078765 


14.0 


1.0010 


0.3231 


0.076952 


14.5 


0.9018 


0.2886 


0.075255 


15.0 


0.8155 


0.2589 


0.073661 


15.5 


0.7398 


0.2331 


0.072161 


16.0 


0.6734 


0.2106 


0.070746 


16.5 


0.6150 


0.1911 


0.069409 


17.0 


0.5629 


0.1738 


0.068143 


17.5 


0.5169 


0.1586 


0.066942 


18.0 


0.4757 


0.1452 


0.065801 


18.5 


0.4386 


0.1331 


0.064715 


19.0 


0.4055 


0.1224 


0.063679 


19.5 


0.3757 


0.1129 


0.062691 


20.0 


0.3488 


0.1043 


0.061747 



TABLE II: Radial velocity and radial (p r ) and tangential (pt) 
components of the black hole momentum as a function of the 
separation in ADMTT coordinates for selected values of the 
separation. The numbers have been produced from a PN- 
inspiral from D = 100M. 



used previously in [ll|. In Table (JTTJ) we tabulate results 
for selected values of the black-hole separation. 



For convenience, we have computed analytical fits for 
the momentum parameters from our numerical solution 
to the PN evolution equations. These fits have been com- 
puted from a PN-inspiral starting at D = lOOAf . With 
the expressions 



2.7069 1.21329 



1.5053 2.60155 



Pr 
Pt 



1.9188 1.76084 



5.3029 9.06417\ 



It 



r 



r 



-2.993 



-3.288 



35.0988 ^ 
± ( fW(r) - r5 . 36702 j 



the relative errors with respect to our numerical results 
are smaller than 0.3% for r and p r and 3.5 x 10~ 4 for 
p t over the range D = 8M to D — 20M. Here P 3 PAr(r) 
is the 3PN-accurate quasi-circular value, which we have 



IV. CONCLUSIONS 

We have presented a conceptually simple method to 
specify very low eccentricity initial-data parameters for 
the numerical evolution of binary systems. We inte- 
grate PN equations of motion from an initial separation 
of D = 40M to the separation we wish to use as the 
starting point for a numerical evolution of the full Ein- 
stein equations. Initial conditions for the PN inspiral 
are taken from a 3PN accurate circular orbit condition. 
These conditions lead to a small initial eccentricity that 
radiates away by the time we read off the parameters for 
our full GR numerical evolution. 

The PN equations are accurate to 3PN order in the 
conservative part when neglecting spins. For spin-spin 
and spin-orbit interactions we have only implemented 
the leading order terms, although these were not used 
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in the application presented here. Radiation-reaction is 
implemented via averaging over orbits, and is accurate 
to 3.5PN order beyond the leading quadrupole contribu- 
tion. The method can in principle be applied to general 
black-hole initial data, and applications to unequal-mass 
and spinning cases will be presented elsewhere. 

We report in detail on an equal-mass inspiral, where 
the evolution of the full Einstein equations is started at 
a coordinate separation of D = 11M (38|. Our method 
reduces the eccentricity by at least a factor of five to a 
value below e = 0.002. Remaining oscillations in the 
coordinate distance of the two black holes cannot clearly 
be identified as eccentricity. We also provide a curve fit 
and table of low-eccentricity parameters for equal-mass 
binaries. 

An important corollary from the success of this method 
is that the input parameters in the initial data construc- 
tion ({mi,Pi, D} in the Bowen-York extrinsic curvature 
and conformal flatness ansatz) actually correspond to the 
physical properties of the black holes during a long-term 
evolution with excellent accuracy. One might instead 
have found that the presence of "junk" radiation spoils 
the initial data, and that after this radiation has left the 
system the dynamics of the black holes are very differ- 
ent from what one would have expected from, for exam- 
ple, a PN evolution with the same physical parameters. 
The fact that PN low-eccentricity parameters translate to 



low-eccentricity numerical evolutions with Bowen-York 
puncture data suggest that neither the junk radiation 
nor the constraint-solution procedure adversely affect the 
physics of the system. 

Pfeiffer, et al. fl"o| have shown that waveforms from 
quasi-circular and low-eccentricity parameters have large 
fitting factors (> 0.99 for I < 4 multipole contributions), 
and conclude that "quasi-circular" waveforms will be suf- 
ficiently useful for gravitational-wave detection. How- 
ever, the use of waveforms with the lowest eccentric- 
ity possible will be necessary to make the most accu- 
rate matching possible to PN inspiral waveforms. Low- 
eccentricity waveforms were compared with PN wave- 
forms in [a], and we will give a further detailed com- 
parison in [39j. 
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